notebook16-tidymodels

Attempting to use tidymodels for processing our IBM krill data and selecting variables of importance for each swimming parameter

Inputs

rf_mod <- ## creates random forest model
  rand_forest(trees = 1000) %>% 
  set_engine("ranger") %>% 
  set_mode("regression")


set.seed(234)  ## comment out to get random runs

rf_fit <- ## fits random forest model to whole dataset
  rf_mod %>% 
  fit(ave.v ~ flow * chl * guano * light, data = df)
rf_fit

rf_pred <- predict(rf_fit, df)  
plot(df$ave.v, rf_pred$.pred, main = "corr - 0.7294325")  ## observed vs predicted
cor <- cor.test(df$ave.v, rf_pred$.pred)  ## gives correlation coef

rf_split <- initial_split(df %>% select(flow, chl, guano, light, ave.v), ## splits cases base on initial results
                            strata = NULL)  ## with strata = NULL splits 11/31, rf_test has all guano absent, 50/50 light split, uneven flow split (0, 3, 3, 5.9, 5.9, 5.9, 8.9 x5), and random chl values.
                                            ## with strata = flow, splits 12/30, rf_test has 1 guano present, 50/50 light split, even flow splits (3x 0, 3, 4x 5.9 and 2x 8.9), and more even chl values.
## play with prop = 0.8, 0.9, etc


rf_train <- training(rf_split)  ##creates training and testing datasets
rf_test  <- testing(rf_split)

## comparisons to test data using ROC and accuracy to measure performance

rf_fit2 <- ## fits random forest model to training dataset
  rf_mod %>% 
  fit(ave.v ~ flow * chl * guano * light, data = rf_train)
rf_fit2

rf_pred2 <- predict(rf_fit2, rf_test) ## compares to test data
plot(rf_test$ave.v, rf_pred2$.pred, main = "corr - 0.2700745")  ## observed vs predicted, can add cor value from cor.test below
cor2 <- cor.test(rf_test$ave.v, rf_pred2$.pred)  ## gives correlation coef
cor2$estimate
### can run with different seeds and splits, strata = NULL
#### loop and run 10 times with diff seeds
## look at consistency of results

Models

library(tidymodels)  # for the parsnip package, along with the rest of tidymodels

── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────── tidymodels 1.2.0 ──
✔ broom        1.0.5      ✔ rsample      1.2.1 
✔ dials        1.2.1      ✔ tibble       3.2.1 
✔ dplyr        1.1.4      ✔ tidyr        1.3.1 
✔ infer        1.0.7      ✔ tune         1.2.0 
✔ modeldata    1.3.0      ✔ workflows    1.1.4 
✔ parsnip      1.2.1      ✔ workflowsets 1.1.0 
✔ purrr        1.0.2      ✔ yardstick    1.3.1 
✔ recipes      1.0.10     
── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────── tidymodels_conflicts() ──
✖ dplyr::combine()       masks randomForest::combine()
✖ purrr::discard()       masks scales::discard()
✖ dplyr::filter()        masks plotly::filter(), stats::filter()
✖ dplyr::lag()           masks stats::lag()
✖ randomForest::margin() masks ggplot2::margin()
✖ recipes::step()        masks stats::step()
• Search for functions across packages at https://www.tidymodels.org/find/

# Helper packages
library(readr)       # for importing data

Attaching package: ‘readr’

The following object is masked from ‘package:yardstick’:

    spec

The following object is masked from ‘package:scales’:

    col_factor

library(broom.mixed) # for converting bayesian models to tidy tibbles
library(dotwhisker)  # for visualizing regression results
library(vip)         # for variable importance plots

Attaching package: ‘vip’

The following object is masked from ‘package:utils’:

    vi

seeing if this fixes error

mean.vel.model <- readRDS(mean velocity.rds)

Error: unexpected symbol in "mean.vel.model <- readRDS(mean velocity.rds"

Node Purity heat map

  ### to save each plot as an individual graph
  #jpeg(filename= paste('~/Bigelow/Figures/tidymodels Output/fit', i,'.jpeg', sep = ''), width = 960,      height = 780)
  #plot(df$response, rf_pred$.pred, main = paste("corr = ", cor$estimate))  ## observed vs predicted
  #dev.off()


  #best_tree <- rf.skill.test(df = df, trees = 2000, "ave.v", prop = 0.5, strata = NULL, do.plot = TRUE) ## not working....
  #bt <- rbind(best_tree, bt)

## look at things outside random forest


rf_training_pred <- 
  predict(rf_fit, rf_train) %>% 
  bind_cols(predict(rf_fit, rf_train, type = "numeric")) %>% 
  # Add the true outcome data back in
  bind_cols(rf_train %>% 
              select(flow, chl, guano, light))

rf_training_pred$light <- as.factor(rf_training_pred$light)
rf_training_pred$guano <- as.factor(rf_training_pred$guano)

rf_training_pred %>%                # training set predictions
  roc_auc(truth = light, .pred...1)

rf_training_pred %>%                # training set predictions, only works for factors
  accuracy(truth = ave.v, .pred...2)

## now that the model has exceptional performance lets move to the test dataset

rf_testing_pred <- 
  predict(rf_fit, rf_test) %>% 
  bind_cols(predict(rf_fit, rf_test, type = "numeric")) %>% 
  bind_cols(rf_test %>% select(flow, chl, guano, light))

rf_testing_pred$light <- as.factor(rf_testing_pred$light)
rf_testing_pred$guano <- as.factor(rf_testing_pred$guano)

rf_testing_pred %>%                   # test set predictions
  roc_auc(truth = light, .pred...1)

rf_testing_pred %>%                   # test set predictions
  accuracy(truth = light, guano)

## differences caused by training set error (bias) by model

###### resampling to the rescue

set.seed(345)
folds <- vfold_cv(rf_train, v = 10)
folds

rf_wf <- ## bundles workflow and random forest model together without a recipe needed
  workflow() %>%
  add_model(rf_mod) %>%
  add_formula(ave.v ~ .)

set.seed(456)
rf_fit_rs <- 
  rf_wf %>% 
  fit_resamples(folds)  ##fits resamples

rf_fit_rs  ## .metrics column contains metrics on model performance

collect_metrics(rf_fit_rs)  ##manually unnests meterics data

rf_testing_pred %>%                   # test set predictions (AS ABOVE)
  roc_auc(truth = light, .pred...1)

rf_testing_pred %>%                   # test set predictions  (AS ABOVE)
  accuracy(truth = light, guano)

Tuning the model

library(glmnet)
library(rpart.plot)  # for visualizing a decision tree
library(vip)         # for variable importance plots

tune_spec <- 
  decision_tree(  ## this is the type of model
    cost_complexity = tune(),
    tree_depth = tune()
  ) %>% 
  set_engine("rpart") %>% 
  set_mode("regression")

tune_spec

tree_grid <- grid_regular(cost_complexity(),
                          tree_depth(),
                          levels = 5)

tree_grid

tree_grid %>% ## shows each level we will tune the model at
  count(tree_depth)

set.seed(234) ## don't understand what these do??
rf_folds <- vfold_cv(rf_train)  ## creates cross-validation folds for tuning

set.seed(345)

tree_wf <- workflow() %>%  ##creates the workflow
  add_model(tune_spec) %>%
  add_formula(ave.v ~ .)

tree_res <- ## resamples and tunes model
  tree_wf %>% 
  tune_grid(
    resamples = rf_folds,
    grid = tree_grid
    )

tree_res  ## gives tuning results

tree_res %>% 
  collect_metrics()  ## collects metrics from tuned models


tree_res %>%
  collect_metrics() %>%
  mutate(tree_depth = factor(tree_depth)) %>%
  ggplot(aes(cost_complexity, mean, color = tree_depth)) +
  geom_line(size = 1.5, alpha = 0.6) +
  geom_point(size = 2) +
  facet_wrap(~ .metric, scales = "free", nrow = 2) +
  scale_x_log10(labels = scales::label_number()) +
  scale_color_viridis_d(option = "plasma", begin = .9, end = 0)

## stubbiest tree with a depth of 1 performed the worst
## deepest tree with depth of 15 did better

tree_res %>%
  show_best(metric = "rmse")  ## shows best model fit

best_tree <- tree_res %>%  ## pulls out data on the best fit
  select_best(metric = "rmse")

best_tree ## summary of best tree model

final_wf <- ## create workflow from best tree model after tuning
  tree_wf %>% 
  finalize_workflow(best_tree)

final_wf

final_fit <- ## create final model from new fit
  final_wf %>%
  last_fit(rf_split) 

final_fit %>%
  collect_metrics()


final_fit %>%  ## plot ROC and compare performance after tuning
  collect_predictions() %>% 
  roc_curve(flow, ave.v) %>% ### NOT WORKING
  autoplot()

final_tree <- extract_workflow(final_fit)  ## extract our final fit for future use
final_tree

final_tree %>%  ## creates workflow plot
  extract_fit_engine() %>%
  rpart.plot(roundint = FALSE)

final_tree %>% ## shows which variables are most important to the model in a plot
  extract_fit_parsnip() %>% 
  vip()

args(decision_tree)

Bigger RF Model

cores <- parallel::detectCores() ## sees how many cores we have to process the data
cores

rf_mod <- ## random forest model generation, parallel processing of models
  rand_forest(mtry = tune(), min_n = tune(), trees = 1000) %>% 
  set_engine("ranger", num.threads = cores) %>% 
  set_mode("regression")

rf_recipe <- ## create random forest model recipe
  recipe(ave.v ~ ., data = df) 
  
rf_workflow <- ## create random forest model workflow
  workflow() %>% 
  add_model(rf_mod) %>% 
  add_recipe(rf_recipe)

rf_mod

val_set <- validation_split(df, 
                            strata = flow, 
                            prop = 0.80)
val_set
# show what will be tuned
extract_parameter_set_dials(rf_mod)

set.seed(345)
rf_res <- 
  rf_workflow %>% 
  tune_grid(val_set,
            grid = 25,
            control = control_grid(save_pred = TRUE),
            metrics = metric_set(rmse))

rf_res %>% ## shows 5 best random forest models out of the 25 candidates
  show_best(metric = "rmse")

autoplot(rf_res)  ## plot results

rf_best <- ## creates model with best predictors
  rf_res %>% 
  select_best(metric = "rmse")
rf_best

rf_res %>% ## collects data for ROC curve plot
  collect_predictions()

rf_best$mtry <- as.integer(rf_best$mtry)


##NOT WORKING
rf_auc <- ## creates set of models with best model and model model for comparison
  rf_res %>% 
  collect_predictions(parameters = rf_best) %>% 
  roc_curve(mtry, .pred) %>% 
  mutate(model = "Random Forest")  

##NOT WORKING
bind_rows(rf_best, rf_res) %>% ## plots model comparisons on ROC curve
  ggplot(aes(x = 1 - specificity, y = sensitivity, col = model)) + 
  geom_path(lwd = 1.5, alpha = 0.8) +
  geom_abline(lty = 3) + 
  coord_equal() + 
  scale_color_viridis_d(option = "plasma", end = .6)

################ last model after tuning 

# the last model
last_rf_mod <- 
  rand_forest(mtry = 8, min_n = 7, trees = 1000) %>% 
  set_engine("ranger", num.threads = cores, importance = "impurity") %>% 
  set_mode("regression")

# the last workflow
last_rf_workflow <- 
  rf_workflow %>% 
  update_model(last_rf_mod)

# the last fit
set.seed(345)
last_rf_fit <- 
  last_rf_workflow %>% 
  last_fit(rf_split)

last_rf_fit

last_rf_fit %>%  ## collect metrics from final model
  collect_metrics()

last_rf_fit %>% ## updates model fit
  extract_fit_parsnip() %>% 
  vip(num_features = 20)

##NOT WORKING
last_rf_fit %>% ## plots best ROC curve, with best set of hyperparameters as predictors
  collect_predictions() %>% 
  roc_curve(ave.v, .pred...1) %>% 
  autoplot()